#####################################################################################
#####			Please cite Georgiadou et al. (2018) as indicated in 			            #####
#####  	https://www.sciencedirect.com/science/article/pii/S026137941830026X		  #####
#####################################################################################

### Note: Set "Replication_Data_Code" as working directory in order to execute code from beginning to end
setwd("~/Dropbox/Research/GRR/Far right comparative paper_2000-2012/submission/PublicationDraft")
rm(list=ls(all=TRUE))

library(ggplot2)
library(ggmap)
library(foreign)
library("rgdal")
library(dplyr) # load dplyr
library(rgeos)
library(broom)
library(maptools)
library(ggpubr)
library(tmap) # load tmap package (see Section IV)
library(viridis)
library(grDevices)
library(RColorBrewer)
library(scales)

setwd("./Data")

farright_data<-read.dta(file="GRR_dataMapsR.dta")

EU_NUTS <- readOGR(dsn = "./NUTS_2010_60M_SH/data", layer = "NUTS_RG_60M_2010")

### Change projection to Google-like Map
proj4string(EU_NUTS)
EU_NUTS <- spTransform(EU_NUTS, CRS("+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs"))
#plot(EU_NUTS,xlim = c(300000,1900100), ylim=c(3500000,11200000))
##EU_NUTS <- spTransform(x, CRS("+init=epsg:4238"))


farright_data$NUTS_ID<-farright_data$nuts_2   #create common ID variable for merging spatial and political datasets
farright_data<-farright_data[which(farright_data$year=="2013"),] #get rid of multiple years (we are plotting averages by nuts 2 through time)

EU_NUTS2013_fortified <- fortify(EU_NUTS, region = "NUTS_ID") ##create fortified databse for plotting in ggplot

farright_data$id<-as.character(farright_data$nuts_2)  ## create character id in political variable for merging with geographical database
map_data <-left_join(EU_NUTS2013_fortified,farright_data, by = c("id", "id"))  ## merge political and geographical databases
map_data$nuts2 <- as.character(map_data$nuts_2)
map_data2 <- map_data[nchar(map_data$nuts2)==4,]
map_data3 <- map_data[nchar(map_data$nuts2)==4,]





##################  Plot ER average 
#Basic ggplot choropleth map for ER
p <- ggplot() +
  # municipality polygons
  geom_polygon(data = map_data2,  aes(fill = er_average, 
                                      x = long, 
                                      y = lat, 
                                      group = group)) +
  labs(x = NULL, 
       y = NULL, 
       title = "", 
       subtitle = "", 
       caption = "")
#Add colours with package "Viridis" and plot only mainland Europe (without French colonies)
p2 <- p+
  scale_fill_viridis(
    na.value = "white",
    option = "magma", 
    direction = -1,
    name="ER vote share",
    guide = guide_colorbar(
      direction = "horizontal",
      barheight = unit(2, units = "mm"),
      barwidth = unit(50, units = "mm"),
      draw.ulim = F,
      title.position = 'top',
      # some shifting around
      title.hjust = 0.5,
      label.hjust = 0.5
    ))+
  theme(legend.position = "bottom",axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),plot.background=element_blank())+borders()+coord_cartesian(xlim = c( -4.5e+06 , 9e+06 ),ylim = c( 3900000 , 12000000 ))
## scale_fill_viridis(option = "magma", direction = -1)
## scale_fill_gradient(low='white', high='black')+





#### Add country borders
map_data3<-map_data[which(nchar(map_data$id)==2),]  ### keep country-level entries in spatial data
p3 <- p2 + geom_path(data=map_data3,aes( 
  x = long, 
  y = lat, 
  group = group))   ### add country borders



##################  Plot PRR average 
#Basic ggplot choropleth map for PRR
q <- ggplot() +
  # municipality polygons
  geom_polygon(data = map_data2,  aes(fill = prr_average, 
                                      x = long, 
                                      y = lat, 
                                      group = group)) +
  labs(x = NULL, 
       y = NULL, 
       title = "", 
       subtitle = "", 
       caption = "")



#Add colours with package "Viridis" and plot only mainland Europe (without French colonies)
q2 <- q+
  scale_fill_viridis(
    na.value = "white",
    option = "magma", 
    direction = -1,
    name="PRR vote share",
    guide = guide_colorbar(
      direction = "horizontal",
      barheight = unit(2, units = "mm"),
      barwidth = unit(50, units = "mm"),
      draw.ulim = F,
      title.position = 'top',
      # some shifting around
      title.hjust = 0.5,
      label.hjust = 0.5
    ))+
  theme(legend.position = "bottom",axis.line=element_blank(),axis.text.x=element_blank(),
        axis.text.y=element_blank(),axis.ticks=element_blank(),
        axis.title.x=element_blank(),
        axis.title.y=element_blank(),
        panel.background=element_blank(),panel.border=element_blank(),panel.grid.major=element_blank(),
        panel.grid.minor=element_blank(),plot.background=element_blank())+borders()+coord_cartesian(xlim = c( -4.5e+06 , 9e+06 ),ylim = c( 3900000 , 12000000 ))
## scale_fill_viridis(option = "magma", direction = -1)
## scale_fill_gradient(low='white', high='black')+

#### Add country borders
map_data3<-map_data[which(nchar(map_data$id)==2),]  ### keep country-level entries in spatial data
q3 <- q2 + geom_path(data=map_data3,aes( 
  x = long, 
  y = lat, 
  group = group))   ### add country borders



setwd("../Graphs")
#### Combine PRR and ER in one plot
Figure1<- ggarrange(q3,p3)
#### Export to pdf
pdf('Figure1.pdf',width = 20,height=8)
Figure1
dev.off()
rm(list=setdiff(ls(), "Figure1"))

setwd("..")






